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Generation of defects and disorder from deeply 
quenching a liquid to form a solid 

A. J. Archer, M. C. Walters, U. Thiele and E. Knobloch 


Abstract We show how deeply quenching a liquid to temperatures where it is lin¬ 
early unstable and the crystal is the equilibrium phase often produces crystalline 
structures with defects and disorder. As the solid phase advances into the liquid 
phase, the modulations in the density distribution created behind the advancing so¬ 
lidification front do not necessarily have a wavelength that is the same as the equi¬ 
librium crystal lattice spacing. This is because in a deep enough quench the front 
propagation is governed by linear processes, but the crystal lattice spacing is deter¬ 
mined by nonlinear terms. The wavelength mismatch can result in significant dis¬ 
order behind the front that may or may not persist in the latter stage dynamics. We 
support these observations by presenting results from dynamical density functional 
theory calculations for simple one- and two-component two-dimensional systems of 
soft core particles. 
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1 Introduction 

Solids with a well-ordered crystalline structure have numerous applications in ma¬ 
terials science. In this article we focus on some of the considerations that determine 
whether a crystal grown from a supercooled liquid will be essentially defect-free 
or will contain substantial disorder that will modify its material properties. We are 
interested in particular in the consequence of quenching a uniform liquid to a tem¬ 
perature at which (a) the thermodynamically stable state is a crystalline solid, and 
(b) the supercooled liquid is linearly unstable with respect to growing density per¬ 
turbations. This occurs when the temperature quench is sufficiently deep. When the 
temperature of the liquid is quenched only a little below the freezing temperature, 
the crystal forms instead via nucleation and growth Kua. However, if the temper¬ 
ature quench is sufficiently deep, then the uniform liquid is unstable and any small 
perturbations in the density grow spontaneously lEHTl. This is often referred to as 
the spinodal regime. One expects that in this case the perturbations will evolve into 
a well-ordered crystalline solid but whether this is the case depends crucially on the 
speed of the crystallisation front that develops from an initial small amplitude den¬ 
sity perturbation Gill. Of course we anticipate that the front speed will be slow for 
shallow quenches and faster for deep quenches, an expectation we confirm here. The 
front speed plays a crucial role since the wavenumber of the density modulations de¬ 
posited behind the crystallisation front is in fact determined dynamically and so may 
differ from the wavenumber of the crystal in thermodynamic equilibrium. We ex¬ 
amine here the physical mechanisms responsible for the speed of the crystallisation 
front and distinguish between the so-called pulled fronts whose speed is determined 
by linear processes and pushed fronts whose speed is determined nonlinearly H. 

The fact that the wavelength of the density modulations deposited by a pulled 
crystallisation front differs in general from the length scale of the equilibrium crystal 
is important in determining the nature of the solid structure that is formed. If the 
length scales are similar, then behind the front few defects are generated as the solid 
relaxes to equilibrium. In contrast, if the resulting wavenumber mismatch is large, 
then the defect density and disorder in the solid are greatly enhanced. We hnd in 
general that this is particularly so for deep quenches where the front speed is high. 

We elucidate these notions by studying the solidihcation of model two-dimensional 
(2D) systems of soft-core particles. We use a simple mean-held dynamical density 
function theory (DDFT) ifToUlTl that has been shown to be quite accurate in repro¬ 
ducing Brownian dynamics computer simulations to demonstrate that a mismatch 
between the wavelength of the modulations deposited behind a solidihcation front 
and the length scale of the equilibrium crystal is indeed likely when the liquid is 
deeply quenched. We also observe that this mismatch results in a disordered struc¬ 
ture behind the front. In systems composed of only one species of particles, the 
disorder subsequently largely disappears as the particles are able to rearrange, heal¬ 
ing the majority of the defects. However, in binary mixtures this disorder may be 
frozen-in. 

In this paper we briehy review the pertinent results from our earlier work Gl^ 
and then focus on characterizing the disorder created by fast crystallisation fronts 
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in several systems of interest and the factors influencing whether the resulting dis¬ 
order persists over long time scales. We emphasise that wavenumber selection via a 
moving front is a linear process only when the temperature quench is deep. Shallow 
quenches to a temperature just below where the liquid and crystalline phases coex¬ 
ist do not generate crystallisation fronts unless the crystal phase is nucleated; the 
resulting front is necessarily a pushed front. Furthermore, this type of fronts persists 
into the linearly unstable state until the quench becomes so deep that the speed of 
the pulled front exceeds the speed of the pushed front |j71. 

Propagation of soft matter crystallisation fronts has been considered before, no¬ 
tably in Refs. Il8l fTSlflhll within the so-called phase-held crystal description in one 
spatial dimension and in Ref. EHISl in two spatial dimensions, in both cases fo¬ 
cusing on the properties of pulled fronts. The more accurate DDFT approach used 
here extends and generalises these results to both types of crystallisation fronts and 
to two spatial dimensions. 

We should mention that the DDFT that we use assumes that the particles in the 
system follow overdamped but stochastic equations of motion - i.e. they are implic¬ 
itly treated as colloidal particles immersed in a solvent that acts as a heat bath. Thus, 
after the quench the system is maintained at the quench temperature. This contact 
with a heat bath eliminates the effects of latent heat release. The absence of particle 
inertia implies that the dynamics are diffusive; no propagating excitations (phonons) 
are possible in such systems. 

This paper is stractured as follows: In ^andl^we describe the origin of the basic 
length scales involved in these processes. In S|4]and|^we summarise and illustrate 
the basic theory determining the speed of pulled fronts. In ^and ^we describe the 
DDFT theory for a single component fluid and the structures generated by crystalli¬ 
sation fronts moving with different speeds, ij^and ^compare and contrast these 
results with the results obtained for binary mixtures. The article concludes with a 
brief summary and a few concluding remarks. 


2 Length scales in liquids and solids 


Consider a system of atoms, molecules, colloids, etc (henceforth referred to as ‘par¬ 
ticles’) which collectively exhibit a phase transition from a liquid to a crystalline 
state. Thermodynamic and structural properties of these two phases can in principle 
be found using classical density functional theory (DFT) lfTjll20H22ll . In DFT it is 
shown that there exists a functional £2[p], together with a minimisation principle 


5f2[p] 

5p 


( 1 ) 


The density prohle p*(r) that solves this equation, i.e. that minimises f2[p], is the 
density distribution of the system at equilibrium. Furthermore, [p*] is the thermo¬ 
dynamic grand potential of the system. Solving Eq. Q for state points in the phase 
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diagram where the liquid is the equilibrium phase yields a density prohle that is 
uniform in space, p(r) = po- In contrast, for the crystal phase, the density prohle 
exhibits a regular array of peaks. From this density prohle, quantities such as the 
crystal lattice spacing i can be determined. 

The functional [p] is highly nonlinear, and quantities characterising the crystal, 
such as £, depend on all terms in the functional. In contrast, quantities such as the 
static structure factor S{k) iflTl of the liquid only depend on the linear response of 
the liquid and so only depend on the terms in [p] that are quadratic in the density 
Huctuation p = p — po. Related to the static structure factor is the Fourier transform 
of the linear density response function ^(r — r'), viz., x{k) = —{po/kBT)S{k), that 
relates the change in the density 5p(r) to a change dVextix) in the external potential 

EQI: 

5p(r) = - j dr'x{r,r')5Vext{.r'). (2) 

This formula applies for both uniform and non-uniform fluids; in particular, in the 
case of a uniform fluid with density po and Vext = 0 perturbed by a small amplitude 
external potential 5Vext{f'\ Eq- <|^ determines the resulting change in the density 
prohle; 5p eee p - po = pQ 

The main point of the above comments is to emphasise that quantities pertaining 
to the crystal, such as £, depend on all terms in [p], but quantities such as S{k) and 
the dispersion relation (a{k) ||2l|8][T2, only depend on the quadratic terms in p. We 
emphasise this point because when a uniform liquid is deeply quenched, the length 
scales of the density modulations that initially grow after the quench are determined 
by (x){k) and so only depend on the quadratic terms in p. In particular, the princi¬ 
pal peak in the dispersion relation (a{k) determines the wavenumber of the density 
Huctuation that grows the fastest. The fact that this wavenumber is determined by 
a quantity that only depends on the quadratic terms in p shows that these fastest 
growing modes need not have the equilibrium wavenumber Injl, i.e. they do not 
necessarily generate the correct density modulations for a perfect equilibrium crys¬ 
tal. From these considerations, one can infer that a deeply quenched liquid may well 
produce a disordered solid, because the length scale of the fastest growing modes 
is not in general equal to i. This argument does not address whether, as solidih- 
cation proceeds, the system can rearrange and subsequently anneal all the defects 
generated in the initial stages of the solidihcation process to produce a perfect crys¬ 
tal. Nonetheless, the observation that the initial dynamics after the quench do not 
in general produce density Huctuations of the correct length scale is an important 
observation. 


* The result in Eq. also applies to non-uniform liquids, i.e. to liquids initially at equilibrium 
with a density profile (r) in an external potential Vom (r) , disturbed by an infinitesimal change 
to the external potential, Vohiit) —> (r). The resulting change in the density profile 5p{r) = 

P»i«v(r) - Pold{r) is then also given by Eq. where 5Vext{r) = V„ew{t) - Void{r)- 
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3 Dispersion relation 

The time evolution of the density distribution p (r, f) in a fluid system of particles is 
given by the continuity equation 


dp 

dt 


= -v j, 


(3) 


where j is the current. This equation is simply a statement of the fact that the number 
of particles in the system is a conserved quantity. To solve this equation of course 
requires an expression for the current j = pu, where u is the local fluid velocity. 
In general, we only have formal expressions for this quantity, and to actually cal¬ 
culate the fluid dynamics, approximations are required. For example, for colloidal 
fluids, where the particles move with stochastic (Brownian) equations of motion, 
the following approximation can be rather good; 


j 


-rpv 


5n[p] 

5p 


(4) 


where F — D/kgT is a mobility coefficient, D is the diffusion coefficient, ks is 
Boltzmann’s constant and T is the temperature. This approximation is the central 
result in DDFT ifTMTl . 

Consider now a fluid with uniform density po with a superposed small amplitude 
density perturbation p = p — po that is either artificially imposed, for example, via 
an external field or a consequence of random thermal fluctuations. Eq. 0 shows 
that the perturbation evolves according to 


dl 

dt 


^p + ff{p^), 


(5) 


where ^ is an operator that is obtained by linearising the full dynamical equation 
([^. For example, in the colloidal case, where Eq. Q is applicable, the operator .if = 
— DpoV^c^^^0, where 0 denotes a convolution, i.e. 0 p = f dr'c^^^(r — 
r')p(r'), andc(^)(r) is the Ornstein-Zernike pair direct correlation function ifT^flTll . 
Linearising Eq. 0 and decomposing p into a sum of different Eourier modes, 

p(r,0=Epke'‘‘"+®W', ^=|k|, (6) 

k 


leads to a dispersion relation for the growth rate co(k) of density fluctuations with 
wavenumber k. In the colloidal case described by Eq. 0 we obtain QISlIIl: 

(o(k) = -Dk^[l-poc(k)J, (7) 

where c(k) is the Eourier transform of Note that for the equilibrium fluid, 

we also have the relation S(k) = [1 —poc(k)]^*. 
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Fig. 1 Typical form of the 
dispersion relation a){k) for 
a uniform liquid that is su¬ 
percooled to the region of 
the phase diagram where it is 
linearly unstable against the 
growth of density modula¬ 
tions with wavenumber k fa q. 



Equation (j^ tells us that wavenumbers k for which 0 }{k) >0 grow while those 
for which (£){k) <0 decay. For a deeply supercooled liquid, it can be the case that 
(o{k) is positive for a band of wavenumbers around the value k — q, say. A sketch 
of the typical form of Co{k) for such a supercooled liquid is displayed in Fig.[^ 
Thus in the early stages following a deep quench, we see the growth of density 
modulations with wavenumber kfa q, leading to the appearance of the length scale 
Inlqin the density distribution of the system. However, as emphasised above, this 
scale may differ from the equilibrium crystal lattice spacing £. 


4 Speed of solidification fronts: marginal stability hypothesis 

A further aspect that has not been mentioned in the discussion so far relates to the 
question of how do solidification fronts propagate into the unstable liquid and what 
length scale density modulations do such solidification fronts produce? In the case 
where the liquid is unstable the properties of the solidification front can often be 
determined from the marginal stability hypothesis ITTI [8l IThl 1^ l24l : Consider the 
leading edge of such a front, where the growing density modulations are still small 
in amplitude and suppose this front is advancing with velocity v. In a reference 
frame that moves with the front, Eq. 0 becomes 

^-l-v-Vp =ifp + ^(p2) (8) 

at 

and so the (complex) growth rate of density perturbations in this moving reference 
frame becomes (Oy{k) = /k ■\ + co{k). The marginal stability hypothesis posits that 
the density modulation at the leading edge of the front has zero growth rate in the 
comoving frame - if this mode were to have a positive growth rate, the front would 
be moving faster than the reference frame; if it were to have a negative growth rate 
the front would retract. In other words, the (complex) group velocity in the comov¬ 
ing reference frame must be zero and the growth rate of the density perturbations 
must also be zero, i.e.. 
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and 


iv + 


d(o{k) 

dk 


= 0 


(9) 


/?e[/k-v + a)(A:)] =0, 


( 10 ) 


where k = kr + iki and v = |v|. These equations can in fact be derived, under appro¬ 
priate conditions, using the method of stationary phase applied to the longtime evo¬ 
lution of an infinitesimal spatially localized initial density perturbation O. Together 
they provide three conditions that are to be solved for the three unknowns, v, kr and 
ki. In the front region, the density profile p{r,t) takes the form Pfront(‘?T)> where 
^ = r — vf , since the moving front generates wavelengths in a periodic fashion. With¬ 
out loss of generality we take the front to move along the x-axis, i.e. v = (v,0,0), 
implying that p(r,t) ~ exp(—k,x)sin(kr(x— vf) + Im[co{k)]t) . Thus kr determines 
the wavelength of the density modulations in the front. More importantly, if no 
phase slips take place, then the wavenumber of the density modulations left behind 
by the front is Il7l 1^ fl^ l24l : 


k* =kr + ^Im[(o{k)]. (11) 

Thus k* is determined by the form of CO{k), which is obtained from the linearised 
equation (j^. Therefore the wavelength 2k/ k* of the density modulation created 
behind the advancing front also differs in general from the equilibrium crystal lattice 
spacing i, just like the wavelength associated with the fastest growing mode. 

Thus even if the length scale determined by the maximum in (£){k) is the same as 
the equilibrium lattice spacing for the crystal, as is the case in the simple PFC the¬ 
ory for the crystal, solidification fronts advancing into a deeply supercooled liquid 
will still generate density modulations with a distinct wavelength, requiring substan¬ 
tial subsequent rearrangements of the system in order to form a defect-free crystal 
without strain. 


5 Marginal stability calculation for simple model 

In this section we perform the marginal stability calculation to obtain the front speed 
V and wavenumber k* of the density modulations created behind the front to show 
how these quantities depend on the dispersion relation. We approximate the disper¬ 
sion relation by making a Taylor expansion around the wavenumber corresponding 
to the principal peak and truncating after the k^ term, 

(0{k) K, (Oi„ —a{k — qf' — b{k — qf’^ ( 12 ) 

where ( 0 „, = (£){k = q) is the maximum growth rate. The coefficient a > 0 is related 
to the width of the principal peak while b measures its asymmetry around the peak 
wavenumber k — q. 
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Substituting Eq. ( [T^ into Eq. (|^ we obtain: 

iv — 2a{d~\-iki) — 2>b{d + ikj)^ = 0, (13) 

where we have written d = k^ — q. Separating the real and imaginary parts of this 
equation, we obtain the following expressions for the front speed and the imaginary 
part of k: 


V = 2{a + 3bd)ki 

(14) 


(15) 


Substituting Eq. ( [T2| into Eq. ( [T0| ) we obtain: 

Re[iv{kr + iki) + ( 0 ,„ — a{d + ikif' — b{d + iki)^] =0 (16) 


giving 

- vki + C0„- a{d^ - kf) - b{d^ - 3dkf) = 0. (17) 

Inserting Eqs. ( [T4) i and ( [T5] l into Eq. ( [T7] i we obtain a cubic equation to be solved for 
d. However, for present illustrative purposes it is instructive to proceed analytically 
on the assumption that d is a small quantity. In this case Eqs. ([T4| and ([T5]l become 


V « 2aki 



(18) 

(19) 


Linearization of Eq. (17 1 in r/ now leads to d = 3b(Om/2a^, i.e. to 

3b(Ofn 




2a^ 


( 20 ) 


This result shows that the wavenumber kr of the density modulation in the advancing 
solidification front is not equal to the wavenumber of the fastest growing mode for 
the quenched uniform fluid, q, unless the peak of the dispersion relation is symmet¬ 
ric, i.e. unless b = Q. We also see that the difference between these two wavenumbers 
grows with O),,,, the magnitude of which is related to the degree of undercooling. The 
deeper the quench, the larger is cu,,,. Moreover, inserting these results into Eq. ([n) 
we obtain the wavenumber k* of the modulations deposited behind the front: 


k* 


b(Of}i 


( 21 ) 


Thus the wavenumber k* differs in general from the fastest growing wavenumber q, 
and neither of these wavenumbers is in general equal to 27r/(' and so defects and dis¬ 
order must be present shortly after a deep quench. Some systems are subsequently 
able to rearrange, but others are not, as we show in the subsequent sections for a 
particular model fluid composed of soft-core particles. 
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6 Model fluid 


We consider a 2D system of soft-core particles interacting via the so-called gener¬ 
alised exponential model of index n (GEM-n) pair-potential; 

= (22) 


where the parameter 0 < e < oo determines the energy penalty for a pair of particles 
to overlap completely, R is the radius of the particles and the exponent n determines 
the ‘softness’ of the potential. When n = 2, the potential varies slowly. In contrast, 
when n is large, as the separation distance r between a pair of particles is decreased, 
the potential increases rapidly from « 0 to a value « £ over a short distance at r = R. 
Here, we consider the case when n = 8. Such soft potentials arise as the effective 
interaction potential between polymers or soft macromolecules in solution ll25l4J7ll . 

Our main reason for considering this model centres on the fact that the structure 
and phase behaviour (i.e. thermodynamics) of this model is well described by a 
rather simple approximation for the free energy. The grand potential of the system 
can be decomposed as follows 111120^23: 

Q[p{r)]=^[p{r)]+ j drp{r){Vext{r)-p), (23) 

where p is the chemical potential and ^[p{r)] = ,^,y[p(r)] + ^exipif)] is the in¬ 
trinsic Helmholtz free energy functional. The ideal gas contribution is 

’^id[p{r)]=kBT J drp{r){ln[p{r)A^]-l), (24) 


where A is the thermal de Broglie wavelength and we use the following mean-field 
approximation for the excess contribution to the free energy 


■^exipir)] = dr'p{r)w{\r-r’\)p{r') 


(25) 


which has been widely used in studies of the structure and phase behaviour of soft¬ 
core systems GlESlISHSll. 

The bulk phase diagrams for the 2D GEM-4 and GEM-8 systems were calculated 
in Ref. |(7| (see also Ref. 1591 ) and are displayed in Fig[^ The diagrams exhibit a 
liquid phase at low densities and/or high temperatures which freezes to form a novel 
cluster crystal phase as the temperature decreases or the density increases. 

We assume that the particles move with overdamped stochastic (Brownian) equa¬ 
tions of motion. Therefore, the time evolution of the non-equilibrium density profile 


p(r,f) may be determined using DDFT, i.e., using Eqs. (23 i - (251 together with 
Eqs. (j^ and (j^. Within the present mean-field approximation, the pair direct corre¬ 
lation function is 
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Fig. 2 The phase diagram of 
2D GEM-4 and GEM-8 fluids. 
The binodals are lines of 
thermodynamic coexistence; 
the spinodals correspond to 
the onset of linear instability 
of the metastable uniform 
liquid. 



cP) 

where j3 = {ksT)^^ 
form; 






yy = -^w(|r-r'|), 


(26) 


5p{r)5p (r' 

and so the dispersion relation (j^ has the following very simple 


a){k) = -Dk^[l +popw{k)], 


(27) 


where w{k) is the 2D Fourier transform of the pair potential w{r). The threshold for 
linear instability of the uniform fluid is determined hy co{k = q) = 0, where q ^0 
is the wavevector at which (o{k) has a maximum, i.e., by 1 + pol5w{q) = 0. This 
leads to a very simple linear density dependence of the onset temperature; ksT — 
\w{k = ^)|po, where the marginally stable wavenumber q at onset is determined by 
the condition 


dw(k) 

dk 



(28) 


For the 2D GEM-8 system the onset wavenumber q Ki 5.26/R. Since v!>(A: = q) ^ 
—0.294-sR^ is independent of the density, the linear instability threshold is a straight 
line in the phase diagram; 

— « 0.294po7;^ (29) 

e 

(Fig. 1^. In addition, the binodals along which the liquid and crystal phases coexist 
in thermodynamic equilibrium also appear to be straight lines in the phase diagram. 
This is not obvious a priori because the binodal calculation requires that one first 
obtains the crystal density profile, which is a highly nonlinear problem. However, fit¬ 
ting the numerically obtained binodals with a straight line proves to be an excellent 
approximation (Fig. |^. For example, for the GEM-8 fluid we And that the binodal 
for the crystal state at coexistence is given by ksT/e « 0.314po/?^ and that of the 
liquid is kgT /£ « 0.339poR^. Thus, when the temperature kgT/e = I, the density 
of the liquid at coexistence is poR^ « 1/0.339 = 2.95 while that of the crystal is 
Pof?2« 1/0.314 = 3.18. 

At thermodynamic coexistence, at the temperature T = (and the chemical 
potential p = Pcoex), a front between the crystal and the liquid state is stationary. 
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However, on decreasing the temperature below T = Tcoex (or increasing ji above 
/r = /icoex)^ the liquid state is no longer the equilibrium state. For a shallow quench, 
the liquid state remains linearly stable but a crystal can still grow if it is nucleated; 
only a crystal seed that is larger than the critical size grows and the interface (front) 
between the two phases advances at a well-dehned speed Vni determined by nonlin¬ 
ear processes 01601. This (pushed) front propagation was studied in detail in Ref. 
01 for the GEM-4 model. At T = Tcoex, the front speed Vni = 0; as the tempera¬ 
ture T decreases below Tioex the speed Vni increases with increasing quench depth 
IT-Tcoexl- 


If the quench is to a temperature T < Tsp, where Tsp is the temperature determined 
by Eq. ( |29| ) at which the uniform liquid becomes linearly unstable (i.e. the spinodal), 
then front propagation via linear processes is possible, with the speed v determined 
by the marginal stability analysis described in (Q However, as can be seen from 
Eqs. ( [T8 | i, ( [T^ and (|20j, v — 0 at T — Tip since — 0 at Tip. ^^s the temperature 
is decreased below Tip, v increases but remains less than Vni for small jT — Tsp]. 
In this regime the front remains a pushed front even though the liquid is already 
unstable 0|6O1|6D. However, the speed v increases faster than Vni with decreasing 
T resulting in a crossover in speeds at temperature T — T^, where Ex < Tp < Tioex- 
At Tx, the two speeds are equal, v = Vni, but for a sufficiently deep quench v > Vni, 
and for these temperatures (T < Ty) it is the linear process that determines how the 
crystal state advances into the unstable liquid 0 [60l [61] . 

The variation of the speed of the crystallisation front with increasing chemical 
potential ji (at fixed temperature) is analogous to that described above for decreasing 
temperature (at hxed chemical potential). The metastable uniform liquid becomes 
linearly unstable at /fsp > /fcoex and for jj. > front propagation via linear pro¬ 
cesses is possible. However, it is only when /r > /ix > Msp that linear processes gov- 




Fig. 3 (a) The front speed v as a function of the chemical potential p for the GEM-4 fluid with 
temperature kgT /e = 1. The solid line is the result of the marginal stability calculation while the 
symbols connected with a dashed line summarize the results from numerical simulations of the 
2D DDFT equations, (b) The wavenumber of the density modulation selected by the moving 
front, the wavenumber k* deposited behind the front (both calculated from the marginal stability 
condition), and the wavenumber k^^ of the equilibrium crystal. The difference between k* and Irgq 
implies that rearrangements behind the front are inevitable as the system seeks to minimise its free 
energy. 
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ern the propagation of the front and the front speed is determined by the marginal 
stability result. In Fig. [^a) (see also Fig. 4 of Ref. [Tl) we show for a GEM-4 fluid 
with temperature hgT /e = \ that for j3/r > j3/rx~21 the front speed obtained from 
solving the DDFT equations numerically in 2D does indeed agree precisely with 
the speed v predicted by the marginal stability analysis. Figure[^b) (see also Fig. 6 
of Ref. 121) compares the wavenumber k* deposited behind a front generated by a 
deep quench (/r > /dx, equivalently T < T^) with the wavenumber corresponding 
to the equilibrium crystal lattice spacing. Our aim below is to explore the conse¬ 
quences of the dramatic difference between k* and keq revealed in the figure for the 
subsequent evolution of the solid phase, and to demonstrate that it is responsible for 
the inevitable presence of defects and disorder. 


7 Structures formed after a quench 


In Fig. 1^ we display a series of snapshots of the density profile calculated for 
a 2D GEM-8 fluid with bulk density PqR^ = 5, quenched to the temperature 
T* = kgTje — I, where the uniform fluid is linearly unstable [cf. Eq. (29 1 ]. The 
times are given in Brownian time units, t* = ksTFR^t, where f* ^ 1 is the time it 
takes for a particle to diffuse a distance ^ R. Solidification is initiated by imposing 
at time f = 0 small amplitude random fluctuations along the line x = 0 on an oth¬ 
erwise uniform density profile with periodic boundary conditions. This produces a 
pair of solidification fronts propagating out to the left and to the right from this line. 
These fronts can be observed in the top left time t* = 1 density profile in Eig. 
Ahead of the front one observes stripe-like oscillations in the liquid, as predicted by 
the marginal stability analysis (see Refs. nziia for other examples of similar fronts). 
Behind this stripe-like precursor the crystal grows but contains defects, a conse¬ 
quence of the mismatch between the dynamically generated length scale 2n/k* and 
the equilibrium length scale £: to lower its free energy and form the equilibrium 
crystal structure, the system must perform significant rearrangements. These rear¬ 
rangements occur locally, generating defects behind the propagating crystallisation 
front. By time t* = 2, the two fronts have collided owing to the periodic boundary 
conditions (see top right panel in Eig.|^ and as time proceeds the crystal density dis¬ 
tribution continues to rearrange, removing defects in order to lower the free energy. 
By time f* = 100 (bottom right) the system is largely occupied by a well-ordered 
hexagonal crystal with the correct equilibrium length scale i. However, even after 
such a long time, there is still a region where the crystal orientation is not aligned 
with the crystal in the centre of the system - i.e. two different crystal grains remain 
present in the system, with defects remaining on the grain-boundary between these 
two. By ‘grain boundaries’ we refer to regions where regular hexagonal ordering is 
absent as occurs along the boundaries of neighbouring regions of misaligned hexag¬ 
onal ordering (the ‘grains’). 

In Eig. 1^ we display a series of density profile snapshots for a system with the 
same density as in Eig.[^ but here quenched to a lower temperature, T* = kgT/e = 
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Fig. 4 The logarithm of the 2D density profile, ln(p(r,r)R^), at times t* = 1, 2, 10 and 100 (going 
from top left to bottom right) obtained from DDFT for PqR^ = 5 and temperature kgT/e = i. 
Plotting the logarithm of the density, rather than the density itself, allows one to see the grain 
boundaries and defects more clearly. The solidification front was initiated at time t* =0 along the 
line X = 0, or equivalently x = 51.2S. 


0.75. In this case the solidification front speed v is faster, so the two fronts have 
already propagated across the system by the time f* = 1. In keeping with the pre¬ 
diction in Eq. (Ill, the mismatch between 2n/k* and the equilibrium length scale ^ 
increases with quench depth (i.e. with the parameter (Om) and so more defects and 
disorder are produced in the system. Over time these mostly anneal out, although 
some persist for very long times. 

To quantify the degree of ordering in the density profiles, we perform a Delauney 
triangulation ||2l[8l[63|63 on the points corresponding to all the maxima in the den¬ 
sity profile. We then calculate the distribution function p(0) of the bond angles 9 in 
these triangles. Perfect hexagonal ordering results in a Delauney triangulation con¬ 
sisting entirely of equilateral triangles, corresponding to a single sharp peak at 60° 
in p(0). In Fig. 1^ we display p{9) at times f* = 2, 10 and 100 after the quench, cor¬ 
responding to the results displayed in Figs.[^and[^ We do indeed see a single peak 
in p{9), centred at 60°. However, the peak also exhibits ‘shoulders’ on either side, 
particularly for short times after the quench, indicating the presence of strain and 
disorder at early times. These shoulders diminish over time, but never completely 
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Fig. 5 Same as 



but for the lower temperature ksT/e 


0.75. 




Fig. 6 The bond angle distribution p{6), calculated from Delauney triangulation on the density 
profiles displayed in Figs. [^and[^ The distribution exhibits a single peak at 60° due to the pre¬ 
dominance of hexagonal ordering although shoulders on either side of this peak are present at 
early times due to the disorder in the system. These decrease over time, as the system rearranges 
to remove defects and lower the free energy. 
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disappear. Comparing the results for the temperature T* = \ (left plot in Fig.|^ with 
those for T* = 0.75 (right plot in Fig. |^, we observe that for shallower quenches 
{T* = 1) the distribution p{d) at time f* = 2 is closer to the later time equilibrium 
distribution than for deeper quenches {T* = 0.75). These results confirm and quan¬ 
tify the impression derived from the 2D density profiles, i.e. that the deeper the 
quench, the greater the disorder created in the system and the more rearrangements 
that are required to equilibrate the system. This is also what one would expect on 
the basis of Eq. GD- 


8 A binary mixture 


The results in the previous section (see also Refs. Ill ID) show that systems com¬ 
posed of only one species of particle, when quenched, are generally able to rearrange 
to form a well-ordered crystalline solid, even if initially there is disorder in the sys¬ 
tem. However, since binary mixtures are far more likely to be glass formers, we now 
consider a binary mixture of soft GEM-8 particles, interacting via the pair potential 
[cf. Eq. @]: 


i{r) = £i_ 


.-{r/R.jr 


(30) 


where the indices i,j= 1,2 label the two different species of particles. We choose 
the pair potential parameters £12/611 = E £22/611 = 1-5, R12/R11 = 1 andR22/7^ii = 
1.3, corresponding to the case where the potentials wii(r) and wnir) are identical, 
but the potential W 22 (r) between pairs of species 2 particles is stronger and of longer 
range. We use DDET to determine the dynamics of the system following a quench 
of the uniform liquid to temperatures where it is linearly unstable and a hexagonal 
crystal is the equilibrium state. We approximate the free energy of the binary mixture 
using a two-component generalisation of Eqs. ( |2^ - (25 1 , input into the DDET for 
mixtures - i.e. the generalisation of Eqs. (|^ and^^ - as described in detail in Ref. 
j. Eor simplicity, we set the mobility coefficients of the two species to be equal. 


i.e.Ti =r 2 = r. 

In Eig.j^we display results for the final equilibrium state obtained from quench¬ 
ing a 50:50 binary fluid with average densities PiRji = P27^ii = 0.2 down to various 
different temperatures, T* = kgT /£\\. The linear instability threshold Tsp for these 
densities is at T* = 0.149. Eor temperatures below this value, the uniform liquid is 
linearly unstable. Above this value, the uniform liquid is linearly stable and for the 
crystal to form, it must be nucleated. Note that at all the temperatures for which we 
display results in Eig.|7] the density profiles of the two species cease evolving after a 
time off* « 100 — 200. In Eig.j^we plot the quantity [pi (r) — p 2 (r)]Rjj and regions 
where pi > P 2 are coloured black while regions where pi < p 2 are coloured white. 
Regions coloured red are where either (i) both densities are small, or (ii) the two 
densities are equal in magnitude, but not necessarily small. It is worth observing 
that around most of the black density peaks there is a pale ‘ring’. This is because 
wherever there is a peak in the density of species 1, there is also normally a peak in 
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Fig. 7 The density difference [pi (r) — p 2 {r)]R[[ after time t* = 100 (when the profiles are almost 
stationary) obtained from DDFT with average densities PiSj j = P 2 ^i i = 0-2 and the pair potential 
parameters ei 2 /fii 1 = l,e 22 /eii = l-5,-Ri2/Sii = 1 andS 22 /Sll = 1.3. The dimensionless temper¬ 
atures T* = kgT /eii are given above each figure. Plotting this quantity reveals the location of the 
species 2 particles (regions coloured white) which aggregate at the defects and grain boundaries. 
The grain size increases with the temperature T*. 


the density of species 2. The peaks in P 2 (r) are, however, broader than the peaks in 
Pi (r) but lower in height, hence the ‘ring’. They are lower in height because many 
of the species 2 particles gather at the grain boundaries; these are the white peaks 
in the profiles in Fig. The eye can easily pick out the grain boundaries in these 
plots, from the lines of white peaks. We also observe that the density peaks at higher 
temperatures are broader than those at the lower temperatures (cf. Lindemann’s cri¬ 
terion for melting ISl). Comparing the different profiles in Fig. [7] we see that the 
deeper the quench, the smaller the size of the crystalline grains and the more disor¬ 
der is present, since there are more grain boundaries in the system. This observation 
is supported by the results in Fig.|^which displays the bond angle distribution p{9) 
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Fig. 8 The distribution p(0) 
of bond angles obtained from 
Delauney triangulation, cor¬ 
responding to the density 
profiles in Fig. [7] Note that 
as the temperature decreases 
p{6) becomes broader, in¬ 
dicating that there is greater 
disorder in the system. 



obtained from Delauney triangulation. For quenches to lower temperatures, p{9) is 
broader, indicating that there is greater disorder in the system. 

These results show that deep quenches to temperatures below the linear insta¬ 
bility threshold result in the formation of states that are highly disordered - a con¬ 
sequence of the fact that the propagating solidihcation front deposits behind it a 
dynamically generated length scale that differs from that of the equilibrium hexag¬ 
onal crystal, thereby introducing defects into the system. In contrast to the one- 
component system discussed in ^ binary mixtures of particles are not able to rear¬ 
range as easily as in a monodisperse system and so the disorder created by the so¬ 
lidihcation front persists for a very long time (indeed indehnitely, within the present 
simple mean-held DFT treatment). 


9 System with competing crystal structures 

We now present results for a different GEM-8 binary mixture, which has pair po¬ 
tential parameters En = £12 = £ 22 . Rn/Rn = 1 and R 22 /R\\ = 1-5. This particular 
system was hrst considered in Ref. IT], and we provide here additional detail. For 
hxed (sufficiently high) total density p = pi + p 2 , where pi and p 2 are the aver¬ 
age densities of species 1 and 2, respectively, the system exhibits several different 
equilibrium crystal structures depending on the concentration <P = pi /p of species 
1. Examples of these are displayed in Eig. All the prohles displayed correspond 
to local minima of the free energy, but we have not checked whether they are the 
global minima at the given state points. At low values of 0 we observe a hexagonal 
crystal, such as that displayed in Eig.j^a). However, for larger values of 0, the sys¬ 
tem can form a binary square crystal structure, such as that displayed in Eig. [^c). 
Alternatively, it can form a binary hexagonal lattice structure, such as that displayed 
in Eigs. I^b) and [^d). Eor high values of 0 the system forms a different simple 
hexagonal lattice, an example of which is displayed in Eig. |^e). In this hexagonal 
crystal the minority species particles occupy the same lattice sites as the majority 
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Fig. 9 Equilibrium crystal structures for the GEM-8 mixture with fiEij = /3e = 1 for all i,j = 
1,2, R 22 /R 11 = 1-5, Rn/Rii = 1 for concentrations (a) <l> = 0, (b) 0.1, (c) 0.25, (d) 0.5, (e) 0.9 
and average total density pSjj = 4. We plot the quantity [pi(r) — p 2 (r)]R[[, so regions where 
Pi (r) > P 2 (r) are coloured black while those where pi (r) < p 2 (r) are coloured white. All profiles 
correspond to local minima of the free energy, but we have not checked whether they correspond 
to global minima at the given state points. In (f) we display the linear stability threshold density as 
a function of concentration <P, for fixed temperature ksT/e — 1. 


species particles, in contrast to the lattice structures in (b)-(d). See Ref. Ill for fur¬ 
ther details. 

We use DDFT to determine the structures that are formed when the uniform 
liquid mixture is quenched. A solidification front is initiated along the line x/Ri\ = 
25 at time f = 0 by adding a small random value along this line. The fact that there 
are several competing crystal structures, in conjunction with the fact that we quench 
to state points that are far from the linear instability threshold (i.e. we perform a 
deep quench) means that the structures that are formed are highly disordered. 

In Fig. [T^ we present results from quenches for three systems with total density 
P^?i=8 but different concentrations: <P — 0.25 (top), <P — 0.5 (middle) and = 
0.75 (bottom)Instead of displaying the density profiles calculated using DDFT, 
we plot the locations of all maxima in the total density profile p (r) = pi (r) + p 2 (r) 
exceeding 50R^^ at the peak. In addition, we colour-code these points, according 
to the local structure. The criteria for this are obtained via a Delauney triangulation 
on the points. Points where the local ordering is square are coloured black, points 
where there is hexagonal order are coloured red and the density peaks with neither 
local ordering are shown as open circles]^ In Fig. 10 we display on the left the 


^ The early time profile displayed in the top right panel in Eig. 12 of Ref. Q shows the corre¬ 
sponding plot for <f> = 0.25 instead of ^> = 0.5 as labelled there. This error is corrected here in 
Fig.[^ 

^ The specific criteria for deciding to which subset a given density peak belongs is as follows: 
After performing the Delauney triangulation on the set of peaks, we consider each triangle. The 








50 

40 

30 

20 

10 

0 

50 

40 

30 

20 

10 

0 

50 

40 

30 

20 

10 

0 

Fig 

X = 

anc 

ai‘e 

the 

doi 

rig] 

are 

an^ 

is r 


of defects and disorder from deeply quenching a liquid to form a solid 


19 



••OO*******0OOOO««OOO*«***0*To#****«OC 
-0 T 0 0000 OOO00* oooo0 00 »0°‘^%00o — 


10 The peaks in the density profile formed by a solidification front initiated along the line 
25 at time f = 0. The system is a GEM-8 mixture with jSe,; = 1 for all i,J =1,2, S 22 /S 11 = 1.5 
^ 12/^11 = 1 and average total density (pi -|-p 2 )R[j = pSn = 8. The results in the top row 
for <P = 0.25, in the middle row for <P = 0.5 and the bottom row for <P = 0.75. In each case 
plot on the left is for an early time t* = 1.6, shortly after the solidification front has exited the 
lain and before the structure has had time to relax, while that on the right is for a later time (top 
t and bottom right are for t* = 40, whilst the middle right is for t* = 400). The density maxima 
;olour-coded according what kind of triangle they belong to in a Delauney triangulation: right- 
led are black, equilateral are red and scalene are open circles. Portions of the hexagonal crystal 
d, whilst the competing square crystal structure is black. 
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Structure formed at time t* = 1 .6, i.e., shortly after the solidification front has exited 
the box, and also the structure at a much later time (at time t* = 40 for <l> = 0.25 
and 0 = 0.75 and the time t* = 400 for the case 0 = 0.5 displayed on the middle 
right in Fig. 10 1 . The structures that are formed are all highly disordered. 

In the 0 = 0.25 case displayed along the top row in Fig. 10 the structure that 
forms initially is a mixture of square and hexagonal ordering, but with a predomi¬ 
nance of squares. There are also regions containing neither structure, largely on the 
grain boundaries. Over time, the regions of square ordering grow at the expense of 
the regions with hexagonal ordering. The density profiles continue to evolve very 
slightly after t* = 40, but the displayed structure is very similar in its statistical 
properties to the final equilibrium structure. 

The bond angle distribution p{9) corresponding to the 0 — 0.25 profiles is dis¬ 
played in the top panel of Fig. [m This distribution has two main peaks at 45° and 
90°, reflecting significant square ordering in the system. Such squares yield right- 
angle triangles in the Delauney triangulation. Initially, at time t* = 1.6, these peaks 
are rather broad, but over time they become sharper, reflecting the presence of grow¬ 
ing domains of a well-ordered square crystal. The hexagonal ordering in the system 
is reflected by the fact that the peak in p{d) at 45° has a ‘shoulder’ on it, extending 
out to 60°. This shoulder diminishes in height over time, but does not completely 
disappear, reflecting the fact that small regions of hexagonal ordering persist. 

In the middle row of Fig. 10 we display the density peaks formed after quenching 
a uniform liquid with concentration 0 — 0.5. Once again the structure that is formed 
contains regions of both square and hexagonal ordering. The size of both of these 
types of domains grows over time. In the middle right of Fig. 10 we display the 
density peaks at time t* = 400, after which the density profiles cease to evolve. The 
corresponding bond angle distribution p{9) is displayed in the middle panel of Fig. 
11 This distribution has three peaks at 45°, 60° and 90°, reflecting the presence 


of a mixture of squares and hexagons. These peaks are rather broad, reflecting the 
significant disorder in the system. 

In the bottom row of Fig. [T^ we display the structure formed from quenching a 
concentration 0 — 0.75 uniform liquid. We see that hexagonal ordering dominates, a 
fact confirmed by the large main peak in p{9) located at 60° displayed in the lower 
panel of Fig. What is particularly remarkable about these 0 — 0.75 results is 
that the peak at 60° is actually sharper at the early times (t* = 1.6) than later times 
(t* = 40): over time the peak broadens! This is due to the fact that in this case 
the solidification front produces modulations with wavenumber k* that is close to 
the wavenumber for the hexagonal crystal structure. However, these do not match 
exactly so that the hexagonal crystal that is initially formed is strained. Over time, 
the system lowers the free energy by introducing defects which alleviate the strain. 
These defects lead in turn to the broadening of the peak in p{9) at 60°. 


comer angles are 6i, 02 and . The triangle is defined as equilateral if 10,' — 0; | <5° for all pairs 
y = 1,2,3. The vertices of these triangles are coloured hlack. Triangles are defined as right-angled 
if for the largest angle 0i we have 10i — 90° | < 5° and for the other two angles 102 — 031 < 5°. The 
vertices of these triangles are coloured red. All remaining vertices which fall into neither of these 
categories are displayed as open circles. 








Generation of defects and disorder from deeply quenching a liquid to form a solid 


21 



0 30 60 90 120 

0 




Fig. 11 The distrib utio n p{6) of bond angles obtained from Delauney triangulation, corresponding 
to the results in Fig.[^ which are for pSj j = 8 and concentration <P = 0.25 (top), rf> = 0.5 (middle) 
and ff> = 0.75 (bottom). In the top plot, the peaks in the distribution at 45° and 90° show that the 
system is largely composed of square ordered local structure. In contrast, the bottom plot, for 
<P = 0.75, has just one main peak in the distribution at 60° showing that the system is largely 
composed of hexagonal ordered local structure. In the middle plot, which is for <P = 0.5, p(0) has 
peaks in the distribution at 45°, 60° and 90° showing that the system contains both square ordered 
local structure and also hexagonal local ordering, as one can also see from Fig.|10| 
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10 Concluding remarks 

In this paper we have discussed a mechanism that results in the formation of disor¬ 
dered structures, when a liquid is deeply quenched to temperatures where the ther¬ 
modynamic equilibrium state is a well-ordered crystal. This occurs because solidi¬ 
fication fronts in deeply quenched liquids propagate via a mechanism that generates 
periodic density modulations in the system with wavelength that is not necessar¬ 
ily the same as the wavelength required for an equilibrium crystal. The wavelength 
mismatch means that the formation of a well-ordered equilibrium crystal state re¬ 
quires significant rearrangements after the front has passed. In monodisperse one- 
component systems, these rearrangements should generally be possible; this is cer¬ 
tainly the case in the model fluid studied here. However, for polydisperse systems 
or multi-component mixtures, such as the binary mixtures studied here, these re¬ 
arrangements are frustrated and in some cases hindered by the fact that there is a 
variety of particle sizes in the system. 

We should emphasise that the front propagation mechanism focused on in this 
paper operates only when the quench is sufficiently deep, to temperatures below 
the crossover temperature T^. Only for T < Tx do solidification fronts propagate via 
the linear mechanism, with the speed v and wavenumber k* determined by linear 
considerations. Above the front speed is determined by nonlinear considerations 
and in this regime the structure formed behind the front is generally much better 
ordered. 

The results presented here are for rather simple 2D model systems composed of 
soft particles. Nevertheless, we believe that the mechanism that we describe should 
be rather general, although much further work is required to determine the nature 
of crystallisation fronts in other deeply quenched systems, and in particular to de¬ 
termine whether one can reach the regime where the fronts propagate via the linear 
mechanism that we describe. In other systems, it may be the case that the speed v 
never overtakes Vni. For example, it is not known whether this regime is physically 
accessible for particles with a hard core; it may be that this regime only arises at 
densities near or even beyond random close-packing. 

As a final point, we should mention that the mechanism described here is not the 
only means of introducing disorder as liquids solidify. Other mechanism include: (i) 
Defects created by impurities, (ii) Different grains, with defects on the grain bound¬ 
aries, generated in the nucleation regime when growing crystals with different ori¬ 
entation nucleated at different points in the system collide, (iii) Defects introduced 
by crystal growth under the influence of external forces or shear, (iv) Disordered 
materials produced in shallow quenches where the growing crystal forms dendritic 
type structures (diffusion limited growth), leading to the formation of crystal grains 
and defects. 
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